﻿using System;
using System.Collections.Generic;
using System.Text;

namespace MLForgeSharp.Models.UnsupervisedLearningModels.Clustering
{
    /// <summary>
    /// 谱聚类
    /// </summary>
    public class SpecClustModel
    {
        private double[][] data; // 数据集
        private int k; // 簇的数量
        private int[] labels; // 每个数据点的簇标签

        public SpecClustModel(int k = 2)
        {
            this.k = k;
        }

        // 计算欧氏距离
        private double EuclideanDistance(double[] a, double[] b)
        {
            double sum = 0.0;
            for (int i = 0; i < a.Length; i++)
            {
                sum += Math.Pow(a[i] - b[i], 2);
            }
            return Math.Sqrt(sum);
        }

        // 构建相似度矩阵
        private double[,] BuildSimilarityMatrix()
        {
            int numSamples = data.Length;
            double[,] similarityMatrix = new double[numSamples, numSamples];

            for (int i = 0; i < numSamples; i++)
            {
                for (int j = 0; j < numSamples; j++)
                {
                    double distance = EuclideanDistance(data[i], data[j]);
                    similarityMatrix[i, j] = Math.Exp(-distance * distance); // 高斯核
                }
            }

            return similarityMatrix;
        }

        // 计算度矩阵
        private double[,] BuildDegreeMatrix(double[,] similarityMatrix)
        {
            int numSamples = similarityMatrix.GetLength(0);
            double[,] degreeMatrix = new double[numSamples, numSamples];

            for (int i = 0; i < numSamples; i++)
            {
                double sum = 0.0;
                for (int j = 0; j < numSamples; j++)
                {
                    sum += similarityMatrix[i, j];
                }
                degreeMatrix[i, i] = sum;
            }

            return degreeMatrix;
        }

        // 计算拉普拉斯矩阵
        private double[,] BuildLaplacianMatrix(double[,] similarityMatrix, double[,] degreeMatrix)
        {
            int numSamples = similarityMatrix.GetLength(0);
            double[,] laplacianMatrix = new double[numSamples, numSamples];

            for (int i = 0; i < numSamples; i++)
            {
                for (int j = 0; j < numSamples; j++)
                {
                    laplacianMatrix[i, j] = degreeMatrix[i, j] - similarityMatrix[i, j];
                }
            }

            return laplacianMatrix;
        }

        // 计算矩阵的特征向量
        private double[,] ComputeEigenvectors(double[,] matrix)
        {
            int numSamples = matrix.GetLength(0);
            double[,] eigenvectors = new double[numSamples, k];

            // 这里使用简单的幂迭代法计算前 k 个特征向量
            Random rand = new Random();
            for (int i = 0; i < k; i++)
            {
                double[] vector = new double[numSamples];
                for (int j = 0; j < numSamples; j++)
                {
                    vector[j] = rand.NextDouble();
                }

                // 幂迭代
                for (int iter = 0; iter < 100; iter++)
                {
                    double[] newVector = new double[numSamples];
                    for (int j = 0; j < numSamples; j++)
                    {
                        for (int l = 0; l < numSamples; l++)
                        {
                            newVector[j] += matrix[j, l] * vector[l];
                        }
                    }

                    // 归一化
                    double norm = 0.0;
                    for (int j = 0; j < numSamples; j++)
                    {
                        norm += newVector[j] * newVector[j];
                    }
                    norm = Math.Sqrt(norm);

                    for (int j = 0; j < numSamples; j++)
                    {
                        vector[j] = newVector[j] / norm;
                    }
                }

                // 保存特征向量
                for (int j = 0; j < numSamples; j++)
                {
                    eigenvectors[j, i] = vector[j];
                }
            }

            return eigenvectors;
        }

        // 训练模型
        public void Train(double[][] data)
        {
            this.data = data;

            // 构建相似度矩阵
            double[,] similarityMatrix = BuildSimilarityMatrix();

            // 计算度矩阵
            double[,] degreeMatrix = BuildDegreeMatrix(similarityMatrix);

            // 计算拉普拉斯矩阵
            double[,] laplacianMatrix = BuildLaplacianMatrix(similarityMatrix, degreeMatrix);

            // 计算拉普拉斯矩阵的特征向量
            double[,] eigenvectors = ComputeEigenvectors(laplacianMatrix);

            // 对特征向量进行 k-Means 聚类
            kMeansModel kMeans = new kMeansModel(k: k);
            double[][] eigenvectorData = new double[data.Length][];
            for (int i = 0; i < data.Length; i++)
            {
                eigenvectorData[i] = new double[k];
                for (int j = 0; j < k; j++)
                {
                    eigenvectorData[i][j] = eigenvectors[i, j];
                }
            }

            kMeans.Train(eigenvectorData);
            labels = kMeans.GetLabels();
        }

        // 获取聚类结果
        public int[] GetLabels()
        {
            return labels;
        }
    }

    // 示例程序
    public class SpecClustModelExample
    {
        public SpecClustModelExample()
        {
            // 示例数据
            double[][] data = new double[][]
        {
            new double[] { 1.0, 2.0 },
            new double[] { 1.1, 2.1 },
            new double[] { 1.2, 2.2 },
            new double[] { 10.0, 10.0 },
            new double[] { 10.1, 10.1 },
            new double[] { 10.2, 10.2 }
        };

            // 创建模型
            SpecClustModel model = new SpecClustModel(k: 2);

            // 训练模型
            model.Train(data);

            // 获取聚类结果
            int[] labels = model.GetLabels();

            // 打印聚类结果
            for (int i = 0; i < labels.Length; i++)
            {
                System.Console.WriteLine("数据点 (" + data[i][0] + ", " + data[i][1] + ") 属于簇 " + labels[i]);
            }
        }
    }
}
